function d = spatial_difference_y(s,dh)
% second order differencing
% indexing: s(j,i)

d = 0*s;

% j = 1, one sided
d(1, :) = ( ...
    -3*s(1,:) ...
    +4*s(2,:) ...
    -s(3,:) ...
    ) / 2 / dh;
% interior
d(2:end-1, :) = ( ...
    s((2:end-1)+1, :) ...
    -2*s(2:end-1, :) ...
    +s((2:end-1)-1, :) ...
    ) / dh^2;
% j = end, one sided
d(end, :) = -( ...
    -3*s(end, :) ...
    +4*s(end-1, :) ...
    -s(end-2, :) ...
    ) / 2 / dh;
